Density functional theory beyond the linear regime: Validating adiabatic LDA 
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We present a local density approximation (LDA) for one-dimensional (ID) systems interacting via 
the soft-Coulomb interaction based on quantum Monte-Carlo calculations. Results for the ground- 
state energies and ionization potentials of finite ID systems show excellent agreement with exact 
calculations, obtained by exploiting the mapping of an iV-electron system in d dimensions, onto a 
single electron in N x d dimensions properly symmetrized by the Young diagrams. We conclude that 
ID LDA is of the same quality as its three-dimensional (3D) counterpart, and we infer conclusions 
about 3D LDA. The linear and non-linear time-dependent responses of ID model systems using 
LDA, exact exchange, and the exact solution are investigated and show very good agreement in 
both cases, except for the well known problem of missing double excitations. Consequently, the 3D 
LDA is expected to be of good quality beyond linear response. In addition, the ID LDA should 
prove useful in modeling the interaction of atoms with strong laser fields, where this specific ID 
model is often used. 

PACS numbers: 31.15.ee, 32.10.Hq, 32.30.Jc 



Over the last years the theoretical description of op- 
tical properties of complex many-electron systems, from 
molecules, to nanostructures and extended systems, has 
re-flourished due to the efficient implementation of time- 
dependent density-functional theory (TDDFT) d Q. 
The good performance shown by the adiabatic local 
density approximation (ALDA) for many finite sys- 
tems has limited the development of exchange-correlation 
(xc) functionals with a more elaborate time-dependence, 
which is clearly in its infancy compared to static DFT. 
However recently many important deficiencies, especially 
of the adiabatic approximation, have been identified [3|- 

a. 

Ultrafast time-resolved optical spectroscopy has re- 
vealed new classes of physical, chemical, and biologi- 
cal reactions, in which directed, deterministic motions 
of atoms have a key role. The advent of free electron 
lasers with attosecond resolution increases the capabili- 
ties of present femtosecond pump-probe experiments, al- 
lowing for a study of the dynamics of non-equilibrium 
electronic systems in real time. In addition, systems of 
all sizes can be investigated, from the atomic scale to the 
most extended molecules (e.g. DNA, proteins and their 
complexes) and solids. Despite those tremendous exper- 
imental advances, the theoretical description of a real 
molecular system subject to ultrashort, intense, and/or 
high-frequency lasers is still in a fledgeling state. Several 
problems need to be addressed, ranging from the non- 



perturbative nature of the physical processes involved 
to the simultaneous description of the (interacting) elec- 
tronic and nuclear degrees of freedom. Therefore, it is 
of paramount relevance to have a theoretical framework 
which allows for a non-perturbative description of elec- 
trons, and at the same time is able to tackle electron- 
ion dynamics in the excited-state. TDDFT seems to be 
the suitable framework to move the realm of density- 
functional methods beyond the linear regime to describe 
the aformcntioncd processes. One important advantage 
is the combined electron and ion dynamics provided by 
TDDFT 0. 

Many physical processes rely on the knowledge of non- 
linear response functions. Therefore, it is very timely to 
provide a systematic study addressing the performance 
of present functionals in the non-linear regime. To assess 
the quality of a functional we need to have appropriate 
data for comparison. Obtaining accurate experimental 
data in the non-linear regime can be very difficult for real 
systems, due to various limitations, e.g. solvent effects 
or additional approximations goin g into the interpreta- 



tion of the collected data [ljj [14J . These problems can 
be avoided by using exactly solvable models which then 
allow for a direct comparison between the exact spetrum 
and an approximate one. Unfortunately, an exact propa- 
gation of even small three-dimensional systems is compu- 
tationally very demanding, and needs further simplifica- 
tion. One possibility is the reduction of dimensionality, 



2 



i.e. the use of one-dimensional (ID) models, where the 
exact diagonalization is feasible as long as the number 
of electrons is sufficiently small. In the present paper we 
work with systems of interacting electrons in ID. Hav- 
ing the exact solution allows us to test orbital dependent 
functionals such as exact exchange (EXX) which can be 
easily transferred to different dimensions. A local den- 
sity approximation (LDA) is achieved, as in the 3D case, 
by quantum Monte-Carlo (QMC) studies of the homo- 
geneous reference, and parametrizing the corresponding 
correlation energy. 

The present work, besides adding fundamental infor- 
mation concerning the relevance of spatial and time non- 
locality in the xc functional, also provides a proper LDA 
paramctrization for electrons interacting via the soft- 
Coulomb interaction in ID systems. This model descrip- 
tion is widely used in the context of high-intensity lasers, 
where above-threshold ionization and high-harmonic gen- 
eration play an important role 15-lij]. Also, ID two- 
electron systems are employed to gain insight into ex- 
act properties of the xc potential and kernel in static 
and time-dependent density functional theory, since these 
systems can easily be solved exactly fl9l - l2~lj ]. 

The ID Hamiltonian for N particles moving in a gen- 
eral external potential v ext reads 
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where Vi n t describes the electron-electron interaction. In 
order to avoid the singularity of the Coulomb interaction 
we employ the soft-Coulomb potential 
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instead. Here, q\ and q 2 describe the charges of the par- 
ticles while a is the usual softening parameter (atomic 
units e = m = % = 1 are used throughout this paper). 
Unless stated explicitly, we use a = 1 for all our calcula- 
tions. Mathematically, it is straightforward to show that 
the Hamiltonian (flj is equivalent to a single particle in 
N dimensions, moving in an external potential consist- 
ing of all the contributions from v ext and v int . The cor- 
responding Schrodinger equation can, hence, be solved 
by any code which is able to treat non-interacting parti- 
cles in the correct number of dimensions in an arbitrary 
external potential. Due to the Hamiltonian being sym- 
metric under particle interchange, Xj O Xh, the solu- 
tions of the Schrodinger equation can be chosen as sym- 
metric or antisymmetric under such an exchange. For 
the simplest case of two interacting electrons both the 
symmetric and antisymmetric solutions are valid, corre- 
sponding to the singlet and triplet spin configurations, 
respectively. However, for more than two electrons one 
needs to separately ensure that the spatial wave function 



is a solution to the TV-electron problem. For example, a 
totally symmetric spatial wave function is a correct so- 
lution for a single particle in TV dimensions, however, for 
TV > 2 there is no corresponding spin function such that 
the total wave function has the required antisymmetry 
to be a solution of the TV fermion problem in ID. We 
solve this problem by symmetrizing the solutions accord- 
ing to all possible fermionic Young diagrams for the given 
particle number TV [22|. The solution of higher dimen- 
sional problems within these symmetry restrictions has 
been implemented into the OCTOPUS computer program 
[23l l2~j ] . Usually, the lowest energy solution is found to 
be purely symmetric and is discarded for TV > 2. With 
increasing number of electrons we also observe an increas- 
ing number of states which do not satisfy the fermionic 
symmetry requirements. 

As a result of reducing the number of dimensions, we 
need to use an appropriate functional for performing the 
DFT calculations. While any orbital functional can eas- 
ily be transferred between dimensions, those function- 
als based on specific systems need to be recalculated. 
This affects the most common functional, i. e. the lo- 
cal density approximation, available only for the normal 



Coulomb interaction in two and three dimensions 25, 2 



an effective Coulomb interaction of a harmonically con- 
fined wire 27, 28|, and some other ad-hoc ID models [29I — 
31[. In this work, we present and use a parametrization 
of the ID LDA obtained from quantum Monte-Carlo sim- 
ulations, which are exact in ID, using the soft-Coulomb 
interaction in Eq. [2j We assess the quality of the approx- 
imation in calculating ground-state properties as well as 
the linear response for various ID model systems. We 
then proceed to calculate the nonlinear response and 
compare the exact one with the ALDA and adiabatic 
exact-exchange (AEXX) spectra. 

The correlation energy of the LDA is parametrized in 
terms of r s and the spin polarization £ = (Nf — TV|J/TV 
in the form 

e c (r s , C) = e c (r s , ( - 0) + C 2 [e c (r s , C - 1) - e c {r s , C = 0)] 

(3) 



with 



£c (r s ,C = 0,l) 



1 r s + Er 2 s 



2A + Br s + Cr\ + Br\ 
xln(l + ar,+/80, 



(4) 



which proved to be very accurate in the paramctriza- 
tion for other ID systems with a different long-range in- 
teraction 27j. Note that the additional factor of 1/2 is 
due the use of Hartree atomic units, as everywhere else 
in the paper. To obtain the exact high-density result, 
known from the random phase approximation [28[, i.e. 
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a 


= 1.0 


a = 0.5 




C = o 


C = i 


C = o 


A 


18.40(29) 


5.24(79) 


7.40(18) 


B 


0.0 


0.0 


1.120(119) 


C 


7.501(39) 


1.568(230) 


1.890(63) 


D 


0.10185(5) 


0.1286(150) 


0.0964(108) 


E 


0.012827(10) 


0.00320(74) 


0.0250(23) 


a 


1.511(24) 


0.0538(82) 


2.431(62) 


P 


0.258(6) 


1.56(1.31) • 10" 5 


0.0142(25) 


m 


4.424(25) 


2.958(99) 


2.922(83) 


av. error 


6.7- 10" 5 


3.3 • 10~ 5 


7.7- 10~ 4 



TABLE I: Values of the LDA correlation energy parametriza- 
tion in Eq. [4] For the most widely used case, i.e. = 1, the 
parametrization is reported for both unpolarized (£ = 0) and 
fully polarized (£ = 1) systems. The error on the last digits 
is given in parenthesis, while the average error (in Hartree) in 
the full density range is given in the last row. 





Exact 


-Etotal 

LDA 


SLDA 


Exact 


IP 

(S)LDA 


(S)LDA 
e HOMO 


H 


-0.67 


-0.60 


-0.65 


0.67 


0.65 


-0.41 


He 


-2.24 


-2.20 




0.75 


0.75 


-0.48 


Li 


-4.21 


-4.16 


-4.18 


0.31 


0.33 


-0.18 


Be 


-6.78 


-6.76 




0.33 


0.35 


-0.16 


He+ 


-1.48 


-1.41 


-1.45 


1.48 


1.45 


-1.18 


Li+ 


-3.90 


-3.85 




1.56 


1.55 


-1.24 


Be+ 


-6.45 


-6.39 


-6.41 


0.83 


0.85 


-0.63 


Li 2 + 


-2.34 


-2.25 


-2.30 


2.34 


2.30 


-2.00 


Be 2 + 


-5.62 


-5.56 




2.41 


2.38 


-2.06 


Be 3 + 


-3.21 


-3.13 


-3.18 


3.21 


3.18 


-2.86 



TABLE II: Total energies and ionization potentials for one- 
dimensional atoms and ions from exact and (spin-)LDA cal- 
culations as well as the eigenvalues of the highest occupied 
Kohn-Sham orbital. All numbers are given in Hartree. 



to leading order in r s , we fix the ratio a/ A to be equal 
to 8/(7r 4 a 2 ) and l/(7r 4 a 2 ) for C = and C = 1, re- 
spectively. In both cases the exponent m is limited to 
values larger than 1. As a result, the number of inde- 
pendent parameters in the function @ is reduced to 7. 
In addition, for a — 1 the denominator can be simpli- 
fied by setting B = 0.0. However, for smaller values of 
the softening parameter a the linear term in the denom- 
inator is important for achieving better agreement with 
the quantum Monte Carlo results. The optimal values 
of the parameters for a = 1 and a = 0.5 are reported in 
Tab.|H and implemented in the OCTOPUS program (23.124 1. 
For more details on the ID QMC methodology and the 



ogy ar 



parametrization procedure we refer to Refs. 

As a first test, we calculated the ground-state energies 
of small atomic systems, for example, a ID helium atom 
with q = 2 in Eq. ^) and two electrons which interact via 
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FIG. 1: Binding energy per atom of the one-dimenstional 
hydrogen molecule as a function of the distance between the 
two ions, exact and LDA calculations for the singlet ground 
state and the first triplet state. 



the soft Coulomb interaction. The ground state energies 
and ionization potentials from the exact and unpolarized 
LDA calculations are given in Tab. HU We include all 
possible systems with one, two, three and four electrons 
in our test. For open-shell systems, we additionally per- 
formed a spin-DFT (SLDA) calculation, where the xc en- 
ergy was spin dependent according to Eq. ([3]) . All atomic 
calculations were performed in a box ranging from -8 to 
8 bohr with a spacing of 0.2 bohr, which ensures the to- 
tal energy to be converged to the accuracy stated in the 
table. 

As we can see, the LDA total energies for the neu- 
tral and positively charged systems agree very well with 
the exact results. As expected, the spin-resolved calcu- 
lations further improve the agreement for the open-shell 
systems. As a result, the ionization potentials, calculated 
as the difference of the total energies of the N and the 
N — 1 electron systems, from the (S)LDA and the exact 
calculations agree almost perfectly. The largest Kohn- 
Sham eigenvalue ^omo onrv partially accounts for the 
total ionization potential, i.e. the ID LDA violates this 
known property of the exact functional [32j. The good 



agreement for the positively charged systems is not re- 
produced for negatively charged ones. For the small sys- 
tems investigated here, LDA does not bind an extra elec- 
tron while the exact calculation shows that the negatively 
charged systems are indeed stable giving total energies of 
-0.73 Ha, -2.35 Ha and -4.17 Ha for H", He", and Li", 
respectively. A comparison with the total energies of the 
neutral systems shows that in the exact calculation the 
additional electron is only very lightly bound in ID. It is 
no surprise that the LDA, with its usual wrong asymp- 
totic behavior of the exchange-correlation potential, does 
not yield stable negatively charged ions. 

As a second test of the new functional we calculate 
the dissociation curve of the ID hydrogen molecule. For 
these calculations we increased the size of the simulation 
box to range from -20 to 20 bohr in order to achieve con- 
vergence also for the stretched molecule. Fig.[T]shows the 
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FIG. 2: Linear (top) and non-linear (bottom) spectra of Be 2+ 
comparing the exact and the ID LDA calculation. The inset 
in the bottom figure shows a zoom into the region from 2.7 
to 3.0 Ha. 



binding energy per atom as a function of the distance be- 
tween the two ions. As one can see, the known pathology 
of 3D LDA is reproduced also in ID. The singlet state 
yields a good description around the equilibrium distance 
of 1.6 bohr with the binding energy being slightly over- 
estimated by LDA. However, the bond breaking is not 
described correctly due to the strong static correlation 
at large distances. The LDA calculation for the triplet 
state yields very good agreement over the whole range of 
distances corroborating the general experience of LDA 
performing better for more polarized systems. 

After having shown that the ID LDA behaves very 
much like its 3D counterpart for ground-state calcula- 
tions, we turn our attention to TDDFT where we use it as 
an adiabatic approximation to the exact time-dependent 
exchange-correlation potential. The propagations were 
performed in a box ranging from -150 to 150 bohr with 
absorbing boundary conditions [23| and a grid spacing of 
0.2 bohr for a total propagation time of 10 3 a.u. 

In Fig.[3]we compare the spectra calculated in a linear 
and non-linear regime from the exact and the LDA cal- 



LUl UJ2 kJ3 IjJ4 IDs LOG UJ7 0.1 O2 O3 

LDA 1.10 1.74 1.90 1.96 2.00 - - 0.22 0.40 - 

EXX 1.13 1.82 2.08 2.20 2.27 2.30 2.32 0.26 0.43 0.52 

exact 1.12 1.81 2.08 2.19 2.26 2.29 2.32 0.28 0.42 0.54 



TABLE III: Excitation energies from linear and non-linear 
response of the ID Be 2+ atom corresponding to the spectra 
in Fig. [2] Excitations from linear response are denoted as u> 
while those from the non-linear spectrum are denoted with O. 
All numbers are given in Hartree. 



culations for a Be 2+ system, i.e. a positive charge with 
q = 4 and 2 interacting electrons in a singlet configura- 
tion. In the linear regime a kick of 10 -4 Ha/bohr was em- 
ployed at t — which was then increased to 0.01 Ha/bohr 
to obtain the non-linear response. The values of the ex- 
citation energies can be found in Tab. IIIII In linear re- 
sponse, we see five peaks in the LDA spectrum which 
compare well with the first five excitations in the exact 
case. As expected, the agreement is better for lower ly- 
ing excitations and gets worse the closer we get to the 
continuum. As a guide for the eye we included the KS 
HOMO energy of the LDA calculation and the exact ion- 
ization potential. The onset of the continuum itself ap- 
pears at too low energies in the LDA calculation missing 
two more clearly visible peaks in the exact spectrum. In 
other words, the LDA fails to reproduce the proper Ryd- 
berg series, a behavior well known from 3D calculations. 
For comparison we also included the results from an EXX 
calculation which shows a slightly better agreement than 
LDA for the first three excitations but, more importantly, 
reproduces the Rydbcrg series due to the correct asymp- 
totic behavior of the corresponding exchange potential. 
The quality of the EXX results also implies that correla- 
tion is of secondary importance in the system for a = 1. 
The non-linear spectrum shows the same excitations as 
the linear spectrum and three additional peaks for the 
exact and the EXX calculation and two additional peaks 
in the LDA spectrum. Their energies are also listed in 
Tab. IIIII Due to the spatial symmetry of the system all 
even order responses are zero and the first non-vanishing 
higher-order response is of third order. The fii = 0.28 Ha 
corresponds to an excitation from the second to the third 
excited state, where the transition from the ground to the 
second excited state is dipole forbidden and, hence, can 
only be reached in a two-photon process. The other two 
frequencies, ^2 = 0.42 Ha and Q3 = 0.54 Ha, correspond 
to the transitions from first to second and second to fifth 
excited state, respectively. Again, both the EXX and the 
LDA calculations yield a good description of the low ly- 
ing excitations, only the third peak cannot be resolved 
in the LDA spectrum. 

One feature of the exact spectrum that is missing from 
both the LDA and the EXX spectra is the small dip at 
2.8 Ha, see inset in Fig. [2] It results from a Fano res- 
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onance |33j, |34j, i.e. the decay of an excited state into 
continuum states. It is missing from both approximate 
spectra due to the double-excitation character of the in- 
volved excited state. Double excitations can only be de- 
scribed in TDDFT if a frequency-dependent xc kernel 
is employed Any adiabatic approximation, however, 
leads to a frequency independent kernel. Hence, double 
excitations, as well as any resulting features, are miss- 
ing from both the ALDA and the AEXX calculations. 
Apart from the well-known shortcomings of not including 
double-excitations and not giving the correct Rydberg se- 
ries, the ID ALDA reproduces both the linear and the 
non-linear exact spectra quite well. 

We have introduced a one-dimensional LDA suitable 
for the description of systems interacting via the com- 
monly used soft-Coulomb interaction. We have shown 
that the one-dimensional functional is of the same qual- 
ity as its three-dimensional counterpart in the calcula- 
tion of ground-state energies of atomic systems and the 
dissociation of small molecules. Also, the linear spectra 
show the same quality known from 3D calculations with 
low energy excitations being well described while Ryd- 
berg and double excitations arc missing. Generally, for 
the ID LDA one can expect the same success and failure 
in applications that are known from 3D calculations, i.e. 
the quality of the LDA results appears to be independent 
of the dimensionality. We emphasize that the ID LDA 
yields a good description not only in linear response but 
also in the non-linear case. Consequently, one can expect 
3D LDA calculations to perform well for the calculation 
of non-linear response, where the experimental data is 
often difficult to interpret. 

The reduced dimensionality of the model systems 
treated in this work allows for a direct solution of the 
interacting problem for small number of particles. The 
comparison between the DFT and an exact calculation 
allows for an assessment of the quality of approximations 
beyond what is possible in three-dimensional systems. 
One-dimensional model systems can provide useful in- 
sight which hopefully will allow for the construction of 
new functionals in the future. 
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